function [P f]=fft1d_pgram(t,y,dopad,dowindow)
%% get periodogram (P) values and the frequency "bins" (f)

% make column vectors
t=t(:);
y=y(:);


%  padding: 0 means no padding, 1 means pad to next power of 2, 
%  2 means pad to nexrt-next power, etc.
if dopad
    % calculate next power of 2 greater than length(t) to pad with zeros
    extra = dopad-1; % add and "extra power" for dopad > 1 (bad style to mix a flag and factor!)
    twopower=2.^((ceil(log(length(t))/log(2))) + extra);
    y = [y; zeros(twopower-length(y),1)];

else % no zero padding
    twopower = length(y);
end

dt=range(t)/(length(t)-1);
dfreq = 1/(dt*twopower);
nyquist = dfreq*(twopower/2); % maxfreq = dfreq*(twopower/2-1);

if dowindow % Optional: window the data with a Hann window 
    w = hann1d((1:twopower)');
    Wss = sqrt(2)*sum(w).^2;
    Y = fftshift(fft(w.*y,twopower)); 
    P = Y.* conj(Y) / Wss; % DFT periodogram - with this normalization, 
                           % the sum over all elements of the spectral
                           % power array equals the variance of the data
else % no windowing    
    Y = fftshift(fft(y,twopower)); 
    P = Y.* conj(Y) / twopower^2; % DFT periodogram - with this normalization, 
                                  % the sum over all elements of the spectral
                                  % power array equals the variance of the data
end

start = ceil(twopower/2)+1;  % ceil() to convert to integer, rounding up (same as length of f, test this!)
P = 2*P(start:end); % Take only first half of spectrum (2nd half is
                           % redundant). Factor of 2 corrects for this.
                           % Last element corresponds to one frequency bin less than the nyquist frequency,
                           % first element to zero frequency (DC) (=data mean)
                         
f = (0:dfreq:(nyquist-dfreq))'; % frequency matrix (bins for Frequ. 0 to 0.5), ' means transposed

P = P(:); % make P a column vector

% Chop off DC component
f=f(2:end); P=P(2:end);

% CH: I moved smoothing into (re)_apply_filters(), which allows single
% or multi-pass filtering and includes the daniell window in addition to
% the rectangular window used here. However, this impacts the signifcance
% lines (see fft1d_theor_signif())

% smooth spectrum (in log space) with a running mean
%P=10.^filter(ones(1,smoothwin)/smoothwin,1,log10(P));

